## this script reproduces the tables and figures (2-4) in the main text ##
## of "Cabinet Durability and Fiscal Discipline" ##
## by David Fortunato and Matt W. Loftis ##

	library(foreign)

## directory ##

	setwd('~/Dropbox/Sovereign bond rates/data/replication')

## first the model summaries ##
## spending models ##
## pooled ##

	pars <- read.dta('spendingPosteriors.dta')
	summary(pars)

	names <- c('Parties in Government', 'ENP', 'Caretaker Time', 'GDP Per Capita',
		'Unemployment Rate', 'Dependency Ratio', 'Trade Openness', 'Spending', 'Maastricht Era',
		'Budgetary Constraint Index (BCI)', 'Parties in Government x BCI', 'Expected Duration', 'Government Ideology', 
		'GDP Per Capita', 'Unemployment Rate', 'Dependency Ratio', 'Trade Openness', 'Intercept')		

	means <- c()
	sds <- c()
	p <- c()
	
	for(i in 1:length(pars[1, ])){
		means[i] <- mean(pars[, i])
		sds[i] <- sd(pars[, i])
		if(means[i] < 0){
			p[i] <- length(pars[, i][pars[, i] > 0]) / length(pars[, i])
		}else{
			p[i] <- length(pars[, i][pars[, i] < 0]) / length(pars[, i])
		}
	}
	
	cbind(names, means, sds, p)

## and fixed ##
	pars <- read.dta('spendingPosteriorsFix.dta')
	summary(pars)

	names <- c('Parties in Government', 'ENP', 'Caretaker Time', 'GDP Per Capita',
		'Unemployment Rate', 'Dependency Ratio', 'Trade Openness', 'Spending', 'Maastricht Era',
		'Budgetary Constraint Index (BCI)', 'Parties in Government x BCI', 'Expected Duration', 'Government Ideology', 
		'GDP Per Capita', 'Unemployment Rate', 'Dependency Ratio', 'Trade Openness', 'Austria', 'Belgium', 'Denmark',
		'Finland', 'France', 'Germany', 'Greece', 'Ireland', 'Italy', 'Luxembourg', 'Netherlands', 'Portugal',
		'Spain', 'Sweden', 'United Kingdom', 'Intercept')		

	means <- c()
	sds <- c()
	p <- c()
	
	for(i in 1:length(pars[1, ])){
		means[i] <- mean(pars[, i])
		sds[i] <- sd(pars[, i])
		if(means[i] < 0){
			p[i] <- length(pars[, i][pars[, i] > 0]) / length(pars[, i])
		}else{
			p[i] <- length(pars[, i][pars[, i] < 0]) / length(pars[, i])
		}
	}
	
	cbind(names, means, sds, p)	

## and now the deficit models ##
## pooled ##

	pars <- read.dta('deficitPosteriors.dta')
	summary(pars)

	names <- c('Parties in Government', 'ENP', 'Caretaker Time', 'GDP Per Capita',
		'Unemployment Rate', 'Dependency Ratio', 'Trade Openness', 'Deficit', 'Maastricht Era',
		'Budgetary Constraint Index (BCI)', 'Parties in Government x BCI', 'Expected Duration', 'Government Ideology', 
		'GDP Per Capita', 'Unemployment Rate', 'Dependency Ratio', 'Trade Openness', 'Intercept')		

	means <- c()
	sds <- c()
	p <- c()
	
	for(i in 1:length(pars[1, ])){
		means[i] <- mean(pars[, i])
		sds[i] <- sd(pars[, i])
		if(means[i] < 0){
			p[i] <- length(pars[, i][pars[, i] > 0]) / length(pars[, i])
		}else{
			p[i] <- length(pars[, i][pars[, i] < 0]) / length(pars[, i])
		}
	}
	
	cbind(names, means, sds, p)

## now with FEs ##
	pars <- read.dta('deficitPosteriorsFix.dta')
	summary(pars)

	names <- c('Parties in Government', 'ENP', 'Caretaker Time', 'GDP Per Capita',
		'Unemployment Rate', 'Dependency Ratio', 'Trade Openness', 'Deficit', 'Maastricht Era',
		'Budgetary Constraint Index (BCI)', 'Parties in Government x BCI', 'Expected Duration', 'Government Ideology', 
		'GDP Per Capita', 'Unemployment Rate', 'Dependency Ratio', 'Trade Openness', 'Austria', 'Belgium', 'Denmark',
		'Finland', 'France', 'Germany', 'Greece', 'Ireland', 'Italy', 'Luxembourg', 'Netherlands', 'Portugal',
		'Spain', 'Sweden', 'Intercept')		
		
	means <- c()
	sds <- c()
	p <- c()
	
	for(i in 1:length(pars[1, ])){
		means[i] <- mean(pars[, i])
		sds[i] <- sd(pars[, i])
		if(means[i] < 0){
			p[i] <- length(pars[, i][pars[, i] > 0]) / length(pars[, i])
		}else{
			p[i] <- length(pars[, i][pars[, i] < 0]) / length(pars[, i])
		}
	}
	
	cbind(names, means, sds, p)

## now, well replicate the figures ##

	betas <- read.dta('spendingPosteriorsFix.dta')

## have a look (duration is #12) ##

	summary(betas)
	quantile(betas[, 12], 0.99)
	dim(betas)

## draw a picture of the effects of going from: ##
## 3 years of expected durability to 1 ##

	effects <- (365 * betas[, 12]) - (3 * 365 * betas[, 12])
	summary(effects)

	pdf('figure2.pdf', width = 10, height = 9)
	plot(1, 1, type = 'n',
		xlim = c(-0.5, 1),
		ylim = c(0, 5),
		xlab = 'Change in Public Spending (%GDP)',
		ylab = 'Density of Effects',
		main = 'Effect of Reducing Expected Duration from 3 Years to 1 Year',
		axes = FALSE)
	axis(1)
	axis(2, at = c(0, 5), labels = c(0, 'Max'))
	support <- c()
	for(i in 1:1000){
		lines(density(effects[(i*100):((i*100) - 99)]), lwd = 0.5, col = rgb(0, 0, 0, .05))
		lines(c(quantile(effects[(i*100):((i*100) - 99)], 0.05), quantile(effects[(i*100):((i*100) - 99)], 0.05)), c(0, 4.5), col = rgb(.6, .6, .6, 0.05))
		if(quantile(effects[(i*100):((i*100) - 99)], 0.05) > 0){
			support[i] <- 1
		}else{
			support[i] <- 0
		}
	}
	lines(density(effects), lwd = 3, col = rgb(0, 0, 0, 1))
	lines(c(quantile(effects, 0.05), quantile(effects, 0.05)), c(0, 4.5), lwd = 3, col = rgb(0, 0, 0, 1))
	text(0.75, 4, paste('Criterion met\n', round(mean(support), 3)))
	text(0.75, 3.5, paste('Overall 2.5 percentile\n', round(quantile(effects, 0.05), 3)))
	effects <- (365 * betas[, 12]) - (3 * 365 * betas[, 12])
	text(0.75, 3, paste('Mean effect = ', round(mean(effects), 3), '% GDP', sep = ''))
dev.off()

## now on to deficits ##
## read in the draws ##

	betas <- read.dta('deficitPosteriorsFix.dta')

## have a look (duration is #12) ##

	summary(betas)
	quantile(betas[, 12], 0.99)
	dim(betas)

## draw a picture of the effects of going from: ##
## 3 years of expected durability to 1 ##

	effects <- (-365 * betas[, 12]) - (365 * betas[, 12])
	summary(effects)

pdf('figure4.pdf', width = 10, height = 9)
	plot(1, 1, type = 'n',
		xlim = c(-0.2, 1),
		ylim = c(0, 4.75),
		xlab = 'Change in Public Spending (%GDP)',
		ylab = 'Density of Effects',
		main = 'Effect of Reducing Expected Duration from 1 Year to -1 Year',
		axes = FALSE)
	axis(1)
	axis(2, at = c(0, 4.5), labels = c(0, 'Max'))
	support <- c()
	for(i in 1:1000){
		lines(density(effects[(i*100):((i*100) - 99)]), lwd = 0.5, col = rgb(0, 0, 0, .05))
		lines(c(quantile(effects[(i*100):((i*100) - 99)], 0.05), quantile(effects[(i*100):((i*100) - 99)], 0.05)), c(0, 4.5), col = rgb(.6, .6, .6, 0.05))
		if(quantile(effects[(i*100):((i*100) - 99)], 0.05) > 0){
			support[i] <- 1
		}else{
			support[i] <- 0
		}
	}
	lines(density(effects), lwd = 3, col = rgb(0, 0, 0, 1))
	text(0.75, 4, paste('Criterion met\n', round(mean(support), 3)))
	text(0.75, 3.5, paste('Overall 2.5 percentile\n', round(quantile(effects, 0.025), 3)))
	text(0.75, 3, paste('Mean effect = ', round(mean(effects), 3), '% GDP', sep = ''))
dev.off()

## and now the the long term spending effects ##
	betas <- read.dta('spendingPosteriorsFix.dta')
	data <- read.dta('replication.dta')
	
	data <- data[data$country == 'Austria' & (data$year >= 1984 & data$year < 2006), ]
	
	medLong <- c()
	medShort <- c()
	hiLong <- c()
	hiShort <- c()
	loLong <- c()
	loShort <- c()
	longTerm <- rep(4:1, 6) * 365
	shortTerm <- rep(2:1, 11) * 365
	
	for(i in 1:21){
		y <- 1984 + i
		if(i == 1){
			long <- betas[, 1] * data$numpar_year[data$year == y-1] + 
				betas[, 2] * data$enp_year[data$year == y-1] + 
				betas[, 3] * data$caretakertime[data$year == y-1] + 
				betas[, 4] * data$gdppercap[data$year == y-1] + 
				betas[, 5] * data$unemp_rate[data$year == y-1] + 
				betas[, 6] * data$dep_ratio[data$year == y-1] + 
				betas[, 7] * data$openc[data$year == y-1] + 
				betas[, 8] * data$outlays[data$year == y-1] + 
				betas[, 9] * data$maasera[data$year == y-1] + 
				betas[, 10] * data$budgetconstraint[data$year == y-1] + 
				betas[, 11] * data$budgetconstraint_inter[data$year == y-1] + 
				betas[, 12] * longTerm[i] + 
				betas[, 13] * data$lr[data$year == y] + 
				betas[, 14] * data$gdppercap[data$year == y] + 
				betas[, 15] * data$unemp_rate[data$year == y] + 
				betas[, 16] * data$dep_ratio[data$year == y] + 
				betas[, 17] * data$openc[data$year == y] + 
				betas[, 33]
			
		
			short <- betas[, 1] * data$numpar_year[data$year == y-1] + 
				betas[, 2] * data$enp_year[data$year == y-1] + 
				betas[, 3] * data$caretakertime[data$year == y-1] + 
				betas[, 4] * data$gdppercap[data$year == y-1] + 
				betas[, 5] * data$unemp_rate[data$year == y-1] + 
				betas[, 6] * data$dep_ratio[data$year == y-1] + 
				betas[, 7] * data$openc[data$year == y-1] + 
				betas[, 8] * data$outlays[data$year == y-1] + 
				betas[, 9] * data$maasera[data$year == y-1] + 
				betas[, 10] * data$budgetconstraint[data$year == y-1] + 
				betas[, 11] * data$budgetconstraint_inter[data$year == y-1] + 
				betas[, 12] * shortTerm[i] + 
				betas[, 13] * data$lr[data$year == y] + 
				betas[, 14] * data$gdppercap[data$year == y] + 
				betas[, 15] * data$unemp_rate[data$year == y] + 
				betas[, 16] * data$dep_ratio[data$year == y] + 
				betas[, 17] * data$openc[data$year == y] + 
				betas[, 33]		
		}else{
			long <- betas[, 1] * data$numpar_year[data$year == y-1] + 
				betas[, 2] * data$enp_year[data$year == y-1] + 
				betas[, 3] * data$caretakertime[data$year == y-1] + 
				betas[, 4] * data$gdppercap[data$year == y-1] + 
				betas[, 5] * data$unemp_rate[data$year == y-1] + 
				betas[, 6] * data$dep_ratio[data$year == y-1] + 
				betas[, 7] * data$openc[data$year == y-1] + 
				betas[, 8] * mean(long) + 
				betas[, 9] * data$maasera[data$year == y-1] + 
				betas[, 10] * data$budgetconstraint[data$year == y-1] + 
				betas[, 11] * data$budgetconstraint_inter[data$year == y-1] + 
				betas[, 12] * longTerm[i] + 
				betas[, 13] * data$lr[data$year == y] + 
				betas[, 14] * data$gdppercap[data$year == y] + 
				betas[, 15] * data$unemp_rate[data$year == y] + 
				betas[, 16] * data$dep_ratio[data$year == y] + 
				betas[, 17] * data$openc[data$year == y] + 
				betas[, 33]
			
		
			short <- betas[, 1] * data$numpar_year[data$year == y-1] + 
				betas[, 2] * data$enp_year[data$year == y-1] + 
				betas[, 3] * data$caretakertime[data$year == y-1] + 
				betas[, 4] * data$gdppercap[data$year == y-1] + 
				betas[, 5] * data$unemp_rate[data$year == y-1] + 
				betas[, 6] * data$dep_ratio[data$year == y-1] + 
				betas[, 7] * data$openc[data$year == y-1] + 
				betas[, 8] * mean(short) + 
				betas[, 9] * data$maasera[data$year == y-1] + 
				betas[, 10] * data$budgetconstraint[data$year == y-1] + 
				betas[, 11] * data$budgetconstraint_inter[data$year == y-1] + 
				betas[, 12] * shortTerm[i] + 
				betas[, 13] * data$lr[data$year == y] + 
				betas[, 14] * data$gdppercap[data$year == y] + 
				betas[, 15] * data$unemp_rate[data$year == y] + 
				betas[, 16] * data$dep_ratio[data$year == y] + 
				betas[, 17] * data$openc[data$year == y] + 
				betas[, 33]		
		}
	
		medLong[i] <- mean(long)
		medShort[i] <- mean(short)
		hiLong[i] <- quantile(long, 0.95)
		hiShort[i] <- quantile(short, 0.95)
		loLong[i] <- quantile(long, 0.05)
		loShort[i] <- quantile(short, 0.05)
	}

## and now make the figure ##
	pdf('figure3.pdf', width = 8)
	plot(1, 1, type = 'n',
		xlim = c(1985, 2005),
		ylim = c(49, 55),
		xlab = 'Year',
		ylab = 'Spending as %GDP',
		main = 'Long-Term Effects of Cabinet Stability',
		axes = FALSE)
	axis(1)
	axis(2)
	points(c(1985:2005), medLong, pch = 17, col = rgb(0, 0, 0, 0.75), cex = 1.5)
	lines(c(1985:2005), medLong, col = rgb(0, 0, 0, 0.75), lty = 3)
	points(c(1985:2005), medShort, pch = 16, col = rgb(0, 0, 0, 0.25), cex = 1.5)
	lines(c(1985:2005), medShort, col = rgb(0, 0, 0, 0.25))
	for(i in 1:21){
		lines(c(1984+i, 1984+i), c(hiLong[i], loLong[i]), col = rgb(0, 0, 0, 0.75), lty = 3, lwd = 2)
	}
	for(i in 1:21){
		lines(c(1984+i, 1984+i), c(hiShort[i], loShort[i]), col = rgb(0, 0, 0, 0.25), lwd = 2)
	}
	c <- seq(1, 21, by = 2)
	text(1995, 49.4, 'Difference in Spending Between Forming 2-year and 4-year Cabinets', cex = 0.9, col = rgb(0, 0, 0, 0.7))	
	text(1984+c, 49, round(medShort[c]-medLong[c], 3), cex = 0.9, col = rgb(0, 0, 0, 0.7))
	legend('topright', legen = c('2-year Cabinet', '4-year Cabinet'), pch = c(16, 17), col = c(rgb(0, 0, 0, 0.25), rgb(0, 0, 0, 0.75)), bty = 'n')
	dev.off()	
